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Abstract 

The present status of extensive air shower (EAS) simulation procedures is 
reviewed. The advantages of combining numerical and Monte Carlo methods 
for the description of EAS development are discussed. Physics content of 
cosmic ray interaction models is briefly described and their predictions are 
compared to the first LHC data. Finally, some outstanding puzzles related 
to cosmic ray composition at the "ankle" energies are analyzed. 



1. Introduction 

Over the last two decades, simulations of extensive air shower (EAS) de- 
velopment have become an important ingredient of experimental analysis of 
high energy cosmic ray (CR) data. The complexity of the corresponding pro- 
cedures is related to the fact that measured EAS characteristics have a very 
indirect relation to the properties of the primary CR particles, resulting from 
a multi-step nuclear-electro- magnetic cascade development. Air shower sim- 
ulations can thus be improved in two directions: i) towards higher precision 
and/or efficiency of the calculations, ii) regarding the physics content. 

Concerning the former, applying EAS simulation procedures to showers 
induced by ultra-high energy cosmic rays (UHECR), one faces the prob- 
lem of enormous calculation time required. Thus, one has to care about an 
optimization of the shower modeling, in order to obtain sufficiently high sim- 
ulation statistics, and about keeping a high accuracy for calculations of both 
average EAS characteristics and their distributions. 

Coming to the latter, of special importance is a correct description of 
the cascade of nuclear interactions of the hadronic component of air showers, 
which acts as a backbone of EAS. On the other hand, the corresponding 
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theoretical models remain just phenomeno logical ones and involve relatively 
large number of adjustable parameters. Thus, further model development 
related both to improvements of the corresponding theoretical description 
and to retuning of model parameters with new accelerator data is desirable. 

2. EAS simulation procedures 

The most transparent approach to EAS simulation is the direct modeling 
of the shower development, as realized e.g. in the CORSIKA program The 
propagation in the atmosphere and interactions of each particle are traced 
using Monte Carlo (MC) methods. The approach has a natural restriction: 
With the number of cascade particles rising proportionally to the energy of 
the primary cosmic ray, so does the computing time required. Hence, exten- 
sion of the procedure to very high energies requires a certain optimization of 
the method, e.g. employing weighted-sampling: Only a number of represen- 
tative particles among all the secondaries produced per interaction is traced 
further by a code; each of those particles acquires thus some weight. 

A classical example is the Hillas's "thinning" method 0], where a sin- 
gle particle per interaction is chosen with the probability E s /E p , E s being 
the energy of the given secondary and E p - the one of the primary particle. 
Correspondingly, the weight of the chosen secondary is w s = E p /E s w p , w p 
being the weight of the primary. Though the method works well for average 
EAS characteristics, it introduces artificial fluctuations in the distributions 
of air shower observables, as discussed e.g. in j^]. The solution to the problem 
was to impose a restriction on maximal weights [3j, in order to reduce the 
magnitude of artificial fluctuations, and to complement the method by an 



"unthinning" procedure [4j. The latter allows one to convert the distribution 
of "weighted" particles coming from the "thinned" EAS simulation to a re- 
alistic particle distribution and to reduce artificial fluctuations to a tolerable 
level. However, the simulation procedure remains time-consuming - as one 
has to find a balance between a sparser "thinning" and not too large weights. 

On the other hand, the efficiency and the accuracy of air shower modeling 
can be significantly improved combining MC and numerical methods, as is 
done in the CONEX Q and SENECA flj codes. Indeed, as EAS fluctuations 
are dominated by interactions and the propagation in the atmosphere of both 
the primary CR particle and of first few generations of the most energetic 
secondaries, one can apply a two-step procedure: Explicit MC simulation 
of the initial stage of the shower and numerical description of secondary 
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hadronic and electro-magnetic (e/m) cascades, based on the solutions of the 
corresponding cascade equationsJll For a number of applications, like calcu- 
lations of fluorescence and Cherenkov radiation profiles of air showers, the 
so-obtained one-dimensional EAS modeling is sufficient. In the more gen- 
eral case, when one is interested in the signal in ground-based detectors, a 
three-step procedure is applied [gj, which includes MC modeling of both the 
highest and the lowest energy part of EAS while intermediate energy range is 
described numerically. Treating most of the cascade with numerical methods, 
one improves the efficiency of the procedure and enhances its accuracy. 



3. Hadronic interaction models 

As discussed in the Introduction, the least certain part of EAS simulation 
procedures is the treatment of hadronic cascades in the atmosphere, which 
involves phenomenological models of hadronic and nuclear interactions. 
Contemporary CR interaction models, like EPOS 0, QGSJET js| and 



QGSJET-II [9j, and SIBYLL lOj, are characterized by a similar physics con 
tent, being designed to treat general hadronic collisions, which involves both 
nonperturbative "soft" and "hard" parton processes. Soft physics is usually 
described within the Reggeon Field Theory framework as soft Pomeron ex- 
changes; hard parton dynamics is treated within the DGLAP formalism and 



implemented in the models either following the so-called miniiet approach [11 
or the qualitatively similar "semihard Pomeron" scheme [12[. Despite these 
general similarities, the models diverge in their predictions, which is both 
due to technical differences in the implementation of the above-discussed ap- 
proaches and, especially, due to different treatments of non-linear interaction 
effects related to parton shadowing and saturation. 

Around the energy of the CR knee, EAS characteristics obtained using 
different models are relatively close to each other. This is both due to the 
model calibration to similar sets of accelerator data and due to ongoing model 



tests by air shower experiments, notably by KASCADE [13j, which resulted 
in serious improvements of certain models. In the discussed energy range, 
the observed EAS characteristics are rather well reproduced by simulations 
under reasonable assumptions on the CR composition, although certain con- 



tradictions persist and composition studies bear model-dependence [14 . 



1 In SENECA, the numerical solution is employed for hadronic cascades while a pre- 
tabulation is used for e/m cascades. 
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The situation changes drastically at higher energies (E > 10 18 eV), where 
model predictions diverge noticeably and certain important observations by 
air shower experiments can not be explained by the present models. 
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4. First measurements at the Large Hadron Collider 

The first LHC data have a strong impact on EAS simulation proce- 
dures and on the interpretation of CR observations. Apart from providing 
additional constraints for hadronic in- 
teraction models, measurements of sec- 
ondary particle production by the CMS 
and ALICE collaborations gave no evi- 
dence on a more rapid energy rise of the 
multiplicity of hadronic collisions than 
predicted by the present CR interaction 
models. This is illustrated in Fig. (TJ 
where the CMS data on the pseudo- 
rapidity density of charged particles in 
pp collisions are compared to the cor- 
responding results of MC models, the 
latter being obtained applying the ex- 
perimental trigger to the MC generated 
hadronic final states. Clearly, the data 
are bracketed by the model predictions. 
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Figure 1: Predictions of CR interaction 
models for dN c h/dri for different cm. en- 
ergies (from up to bottom: y/s — 7, 2.38, 
0.9 TeV) compared to CMS data 0. 



5. UHECR puzzles 

Historically, the first strong indication on a much higher EAS muon con- 
tent than predicted by simulations has been observed in the HiRes-Mia anal- 
ysis [16(] for CR energies Eq > 10 17 eV. More recently, the Pierre Auger 
collaboration has demonstrated using 4 different methods that experimen- 
tal data favor a much higher (by a factor ~ 1.5) number of muons at 
ground than predicted e.g. by the QGSJET-II model [171 ] . Such a strong 
enhancement can not be achieved with the present CR interaction models, 
as it would require to increase the multiplicity of proton-air and pion-air 
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collisions by an order of magnitude over a wide range of energies [1 
On the other hand, the data are marginally consistent with simulation re- 
sults for iron-induced air showers [17]. As the measurements refer mainly to 
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primary energies below 10 EeV, the required change to an iron-dominated 
composition has to occur around the ankle of the CR spectrum. 

Another striking result obtained was the sharp decrease of the width of 
the shower maximum distribution RMS(X max ) at E > 10 18 eV [20j. As dis- 



cussed in 2lj], RMS(X max ) is an almost model-independent measure of the 
CR composition. Indeed, for proton-induced EAS this quantity is mainly de- 
fined by the mean free pass of the proton X p ~ l/cj£!ii r , which sets the lower 
limit on the RMS(X^ ax ) around 50 g/cm 2 . Fluctuations related to the geom- 
etry of p-air interactions (higher /smaller inelasticity for "central" /peripheral 
collisions) can only increase the corresponding value. On the other hand, 
in case of Fe- induced EAS, RMS(X^ X ) is rather dominated by the fluctua- 
tions of the collision geometry, primarily, via the variations of the number of 
"wounded" projectile nucleons (which participate in particle production) and 
via the fragmentation of the nuclear spectator part [22]. Even for extreme 
assumptions, RMS(X 1 ^! lx ) can not exceed some 30 g/cm 2 . 

The observed decrease of RMS(X max ) from ~55 g/cm 2 at 1 EeV to ~30 
g/cm 2 at 30 EeV may be naively interpreted as a change from a p-dominated 
to an Fe-dominated com pos ition in the discussed energy range 
as correctly noticed in 
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20j . In reality, 



05 



adding 

a small admixture of iron nuclei to 
the pure proton composition increases 
the width of the distribution. As 
an illustration, in Fig. |2] we compare 
the Pierre Auger data on RMS(X max ) 
with EAS simulation results (using 
QGSJET-II) considering a simple 2- 
component (p + Fe) composition and 
assuming the partial abundances j\ to 
change smoothly between 1 and 30 EeV: 
f p (E) = f p (l EeV) [1 - lg(£/EeV)/1.5], 
f Fe (E) = 1 - f p (E). Comparing the 
cases / p (l EeV) = 1 and f p {l EeV) = 
0.4, we see that it is the latter choice 
which is supported by the data, i.e. iron 
dominance at the ankle is favored. 

On the other hand, a heavy CR composition in the EeV energy range is 
at variance with HiRes and Pierre Auger data on the average X" max , which 
are consistent with the proton-dominance. In order to reconcile the latter 




Figure 2: RMS(V max ) for an energy- 
dependent CR composition (as described 
in the text): / p (l EeV) = 1 (solid) and 
/ p (l EeV) = 0.4 (dashed); points - 
Pierre Auger data [201 ] . 



5 



with a heavy composition, a much deeper (than presently predicted) shower 
penetration in the atmosphere has to be assumed. In practical terms this 
would require a factor of 2 decrease for cr^f^r compared to the current model 
predictions [19]], which would correspond to a similar reduction for the total 
pp cross section at y/s ~ 6 TeV, i.e. only slightly above the Tevatron energy. 
Such a sudden fall down of 0"*°* would imply very exotic physics. 

6. Conclusions 

Significant progress in the modeling of air showers has been achieved 
in recent years. Combing numerical and MC methods, as realized in the 
CONEX and SENECA codes, allowed one to increase the efficiency and the 
accuracy of EAS simulations. The description of hadronic collisions by the 
corresponding MC models has been considerably improved, particularly, con- 
cerning the treatment of non-linear interaction effects. Predictions of various 
hadronic MC models have converged to each other, both due to the improved 
theoretical description and due to model tests by EAS experiments. 

The first LHC data on secondary hadron production in pp collisions pro- 
vide no indication on a more rapid energy rise of the multiplicity than pre- 
dicted by the present CR interaction models. Thus, a strong rise of EAS 
muon content at 1 4- 100 PeV is not supported by the collider observations. 

The strong enhancement of EAS muon content and the sharp decrease 
of the width of A max distributions, as observed by the Pierre Auger collab- 
oration above 1 EeV, can not be explained in the framework of the present 
CR interaction models, unless an iron dominance of the CR composition at 
the ankle is assumed. The latter assumption is however in a strong conflict 
with X max measurements at 1 -j- 10 EeV by all the experiments working in 
this energy range. In turn, the data on average A max can not be reconciled 
with a heavy CR composition without invoking very exotic physics. 
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